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The order of the superconducting phase transition is a classical problem. Single-component type-2 
superconductors exhibit a continuous “inverted-XY” phase transition, as was first demonstrated for 
17(1) lattice London superconductors by a celebrated duality mapping with subsequent backing by 
numerical simulations. Here we study this problem in multiband 17(1) London superconductors and 
find evidence that by contrast the model has a tricritical point. The superconducting phase transition 
becomes first-order when the Josephson length is sufficiently large compared to the magnetic held 
penetration length. We present evidence that the huctuation-induced dipolar interaction between 
vortex loops makes the phase transition discontinuous. We discuss that this mechanism is also 
relevant for the phase transitions in multicomponent gauge theories with higher broken symmetry. 


The problem of the order of the superconducting phase 
transition beyond mean-field approximation is more than 
40 years old. In [T] it was observed that if phase fluctua¬ 
tions are neglected in the Ginzburg-Landau model, fluc¬ 
tuations of the vector potential make the superconduct¬ 
ing phase transition first-order (see also [2]). This sce¬ 
nario applies only for strongly type-1 superconductors, as 
was demonstrated by later works [M5]. These works con¬ 
sidered a lattice London model, neglecting fluctuations 
of the modulus of the order parameter but taking into 
account phase fluctuations. It was demonstrated that 
the superconducting phase transition beyond mean-field 
approximation is caused by the proliferation of vortex 
loops with short-range interaction set by the magnetic 
field penetration length A. For low temperatures only 
small vortex rings are thermally excited, while at the 
thermally excited vortex loops loose their line tension 
and extend throughout the entire system. 

In Refs. [3H5] a duality mapping was established be¬ 
tween the statistical problem of the normal state above 
Tc in a superconductor and a statistical description of 
a superfluid state below such that the temperature 
axis is reversed, see also [S]. Thus the phase transi¬ 
tion in the lattice London superconductor is called the 
“inverted-3D XY” transition. The London limit approx¬ 
imation is justifiable for extreme type-2 superconductors. 
The question at what parameter values of the Ginzburg- 
Landau model the inverted-3D XY phase transition turns 
into first-order has subsequently been studied using var¬ 
ious approaches. Some of the attempted analytical ap¬ 
proaches [7] suggested that the continuous phase transi¬ 
tion extends slightly into the region where the Ginzburg- 
Landau parameter k is smaller than I (or, equivalently, 
in traditional units smaller than 11^/2). However this 
approach was based on uncontrollable assumptions, and 
thus numerical simulations were required to address that 
question. Early numerical work were consistent with 
the existence of a tricritical point [8], as well as one- 
loop renormalization-group calculations [5] . The largest- 


scale Monte Garlo simulation performed so far [TO] sug¬ 
gests that the continuous phase transition indeed extends 
slightly into the region of k < I. All these works sug¬ 
gested that in London model (which corresponds to tak¬ 
ing extremely type-2 limit) the phase transition is con¬ 
tinuous. 

Recently there has been a surge of interest to mul¬ 
ticomponent 17(1) X C7(I) and SU{2) generalizations of 
this problem in the context of multicomponent super¬ 
conducting condensates HM and of deconfined quan¬ 
tum criticality proposals where such Ginzburg-Landau 
models were argued to arise as an effective held theory 
[HHiii. Initially it was suggested that SU{2) as well 
as 1/(1) X t/(l) superconductors with equal phase stiff¬ 
nesses possess a continuous phase transition from a state 
with fully broken symmetries to a normal state in a novel 
universality class [nmgiis], but later analysis revealed 
hrst-order phase transitions [m na uMi] in such sys¬ 
tems. The flowgram method proposed in |21j shows that 
the phase transition remains hrst order in the limit of 
inhnitesimally small electric charge. First-order phase 
transitions were also found in L7(l) x Z 2 systems, where 
a Z 2 symmetry associated with broken time-reversal sym¬ 
metry is caused by a frustrated intercomponent coupling 

m- 

Most of the superconductors which are of great current 
interest are multiband superconductors [2HU2H] . In the 
London limit multiband superconductors are described 
by several phases coupled by a Josephson interaction [55]. 
Here we consider the problem of the phase transition in 
a multi-component London superconductor with conden¬ 
sates tfa = IV’al exp(i0(“)) (a = 1,2,...), with constant 
\ipa\, described by the free-energy density 

/ E + eA)2 + ^ (V X A)2 

a 

-^Vab\'lpa\\'lpb\cOs{e^°-'> -0^^^'^), ( 1 ) 

a>b 

where e is the electric charge coupling and rjab determines 
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the strength of the Josephson phase-difference-locking 
term between bands a and b. We are interested here in 
the case where the Josephson coupling breaks the sym¬ 
metry explicitly to 17(1), and to restrict the parameter 
space we will consider rjab = rj- EH] a duality map¬ 
ping was discussed for a class of 17(1) multi-band mod¬ 
els. Namely, it was discussed that the dual version of 
the model is described by proliferation of a single kind of 
directed loops. That indicates that if the phase transi¬ 
tion is continuous then it should be in the “inverted-3D 
XY” universality class. However it was established re¬ 
cently that intermediate-length-scale interaction in the 
directed-loops model can make the phase transition first- 
order [30] . Thus the origin of the phase transition in 
multi-band models requires a careful study. 

We begin by examining the two-band case, then with¬ 
out loss of generality rj can be taken to be non-negative, 
such that in the ground state 6*7) — 6*7) = 0. From 
the free energy two characteristic length scales can be 
identified: the first one is the London magnetic penetra¬ 
tion depth A and the second one is the Josephson length 
the latter at which the system restores the ground 
state from small deviations in the phase difference. As 
previously mentioned, the London model is justified for 
strongly type -2 multiband superconductors, where the 
coherence lengths associated with the densities are much 
smaller than A and ^j. The Josephson length can be 
identified by writing 6*71(2;) _ 6*71(2:) = S{x), imposing 
the condition J(0) = So, and expanding the Josephson 
term. Then the system recovers ground state value of the 
phase difference according to the exponential law 6{x) = 
(5oexp(-a:/0) with = ^/|V’i||^ 2 |/(f?(IV'iP + l^ 2 p))- 
The standard expression for the magnetic field penetra¬ 
tion depth is A = 1/ 

When 77 = 0, the model has 17(1) x C7(l) symmetry. 
As mentioned above, the phase diagram of such a sys¬ 
tem have been previously studied, in two dimensions 
three dimensions in an external field [H EH [31] and 
three-dimensional cases without an external field mm]. 
In three dimensions, at a low but non-zero value of the 
coupling constant e, the model exhibits a single first- 
order phase transition from the C/(l) x 17(l)-state to the 
normal state. For large electric charges there are two 
phase transitions. The reason for occurrence of the sec¬ 
ond phase transition is the following: at a lower critical 
temperature a proliferation of bound states of vortices 
with winding in both phases takes place: i.e. vortices 
for which A 07 ) = V07) = 2tt, A 07 ) = V07) = 27r 

where the integration path cr encloses two cores, such vor¬ 
tices are denoted by (1,1). The “elementary” vortices in 
the individual condensates are bound into these compos¬ 
ite objects through the coupling to the vector potential. 
In the resulting state the individual phases 6*71 a^e disor¬ 
dered but the phase difference is ordered mmEiiEi]. 
This state is called a metallic superfluid, paired state or 
super-counter-fluid. At elevated temperatures the tran¬ 


sition into the normal state is driven by the proliferation 
of individual (also termed fractional) vortices A6*7) = 
£^V 07 ) = 27r, A 07 ) = £,.V6*7) = 0 , denoted by ( 1 , 0 ) 
and A6»7) = £ V6»7) = 0, A6*7) = £ V6»7) = 27r, de¬ 
noted by (0,1). For the properties of the individual and 
composite vortices in this model see mm- 

It has been discussed that the existence of the paired 
states in multicomponent gauge theories could be related 
to first-order phase transitions m]- Indeed, even at the 
level of mean-field analysis the direct 17(1) x 17(1) or 
517(2) phase transitions are first-order in the vicinity 
(on the phase diagram) of the paired state. However 
mean-field analysis is qualitatively wrong for this prob¬ 
lem in general, failing to capture first order phase tran¬ 
sition happening in the e —^ 0-1- limit. By contrast, the 
model Q shares the same symmetry as single-component 
superconductors, indeed the Josephson coupling breaks 
the symmetry explicitly to 17(1). 

By discretizing the model Q onto a cubic grid of size 
L, rewriting the gradients with finite differences, the ki¬ 
netic energy terms into XY-model cosines, and inserting 
the previously identified length scales, we obtain the lat¬ 
tice Hamiltonian for the two-band case 


i 


E E cos(6»,^“l - 6*1“^ -t 

a fi 


1 iv^inv^2p 
ej w + w 


cos(0|^^ 



( 2 ) 


where i is a lattice site index, /i = x,y,z and [A x = 
Ai^' J- Ai^rr is the discrete lattice 

curl, where we have denoted x' = y, y' = z and z' = x. 
We simulate the temperature-driven phase transitions of 
the Hamiltonian ([^ with a parallel tempering Metropolis 
Monte Carlo algorithm |33H53] and use periodic bound¬ 
ary conditions. The trial moves consist of selecting a site 
at random and proposing new values for all five site de¬ 
grees of freedom (two phases and three vector potential 
components). The phases are updated by drawing ran¬ 
dom numbers between 0 and 27r, and the components of 
the vector potential are updated by adding (either posi¬ 
tive or negative) random numbers to the old values. Par¬ 
allel tempering swap trial moves between adjacent tem¬ 
peratures are attempted with a fixed frequency jSS] [37] of 
once every sweep. To collect data we perform simulations 
of at least 10 ® sweeps and we perform the simulations 
with linearly distributed temperatures. 

We begin by determining how a finite Josephson cou¬ 
pling term affects the maximal value of the heat capac¬ 
ity and the corresponding temperature. We consider the 
case of twin bands = \^ 2 \'^ = 1 with A^ = 0.2 

and vary 77 over several orders of magnitude. Fig. 
shows the maximal value of the calculated heat capacity 
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FIG. 1. The maximum of the calculated heat capacity (left 
vertical axis, blue) and its corresponding critical temperature 
(right vertical axis, red) for L = 24 (circles) and L = 36 
(triangles), with Josephson conpling rj varying over several 
orders of magnitnde. 


as well as the corresponding temperature, taken to be 
the critical temperature for two system sizes. Each point 
on the curves corresponds to a simulation and the lines 
are guides to the eye (note however that the tempera¬ 
ture curves overlap). The simulations have a fixed set 
of temperatures each, chosen manually to cover the heat 
capacity peak. 

As is seen in Fig. the heat capacity peak is larger 
in the region where the Josephson interaction is small. 
As the Josephson interaction is increased, the Joseph¬ 
son length becomes smaller as does the magnitude of the 
heat capacity peak. In the right-most region where the 
Josephson length is smaller than the magnetic penetra¬ 
tion depth, the magnitude of heat capacity peak is not 
visibly different for the two system sizes. The critical 
temperature raises as the heat capacity peak diminishes 
and is not visibly affected by changing the system size. 
The behavior of the displayed quantities suggests that 
the left-most and right-most points of Fig. represents 
two distinct physical regimes. 

For the left-most point in Fig. (with rj = 0) the 
symmetry of the system is 17(1) x 17(1), which, as pre¬ 
viously mentioned, is known to have a first-order phase 
transition [191 HIl [31] • This, combined with the observa¬ 
tion that the heat capacity maximum grows substantially 
with the system size in the region of finite and small rj, 
suggests that the phase transition is of first-order also for 
finite rj in the region with X < ^j. In the region of A > 
however, any growth in the heat capacity maximum with 
system size is not visible in Fig. |l] suggesting there is 
a tricritical point and a continuous phase transition for 

^ > O- 

We consider first the scaling properties for r] = 10“^ 
(with « 23) which is in the A < ^j-region of Fig. 
The internal energy histogram at the phase transition is 
seen in Fig. a) to be bimodal, and the bimodality is 
enhanced by increasing the system size. The free-energy 
barrier AF = ^n{Pmax/Pmin) {/3c is the inverse 
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FIG. 2. For 77 = 10“^ the scaling behavior of the energy 
histograms (normalized distributions of the internal energy 
per site U/L^) at the phase transition is indicative of a first- 
order phase transition. The histograms for various lattice 
sizes are shown in a), the free-energy barrier AF is shown in 
b) and the latent heat AH in c). 


critical temperature and P the energy distribution) is 
seen in Fig. [^b) to be proportional to and the latent 
heat AH (i.e. the distance between the peaks) is seen in 
c) to not vanish, indicating a first-order phase transition 
EiiD]. By contrast, we observed no bimodal features 
or tendencies for any system with A > ^j. 

We find that when the Josephson interaction is in¬ 
creased, the system size required for a bimodal energy 
distribution to develop also increases. This makes it com¬ 
putationally difficult to reproduce the scaling results of 
Fig. [2] for systems where the Josephson length is much 
smaller than the box length. To further test these phe¬ 
nomena we consider a three- and four-band models with 
\4’i\’^ = 1 with a Josephson coupling which locks all the 
phase differences to zero in the ground state. In such 
a system we find that the first-order phase transition 
becomes stronger and the measurement of bimodal en¬ 
ergy distributions for ^ L is possible. In Fig. is 
shown the first-order scaling of the energy histograms for 
a three-band model with 77 = 0.03 such that « 3.3 up 
to L = 44, more that ten Josephson lengths. The results 
of Fig. [^are 10® sweep simulations with A^ = 0.1. In 
Fig.i is shown simulation results for a four-band model 
with an even stronger Josephson coupling of 77 = 0.07, 
and A^ = 0.075. The four-band results are 2 • 10® sweep 
simulations with sampling during the last 500000 sweeps. 
With the growing number of components we find stronger 
signatures of the first order phase transitions. 

Consider next the scaling properties for the two-band 
model at the point 77 = 10 ® ® which is in the region of 
A > in Fig. (with « 0.40). In Fig. |^it is seen 
that the energy cumulant Vl = 1— (U^) / {3{lP)‘^) has for 
77 = 10 “® a minimum that is persistent against scaling, 
as it should for a first-order transition [IT]. For 77 = 10® ® 
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FIG. 3. First-order scaling of the energy histograms for the 
three-band model also hold when <C L, here = 3.3 with 
scaling up to L = 44. 
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FIG. 4. First-order scaling of the energy histograms is more 
pronounced for a four-band model. 


however, the energy cumulant has no distinct minimum, 
a behavior which is consistent with a continuous phase 
transition where the energy cumulant reaches the triv¬ 
ial limit of 2/3 everywhere in the thermodynamic limit. 
Indeed, in the limit of strong interband coupling the sys¬ 
tem should recover the standard 3D inverted-XY phase 
transition like in a single-component model HE]- 

We interpret the results as suggestive to an emer¬ 
gent attractive interaction between vortex lines. It has 
been recently demonstrated that in a single-component 
directed-loops (j-current) model, a modification of the 
vortex interaction potential from short-range repulsive 
to short-range non-monotonic asymptotically attractive, 
causes a conversion of the inverted-XY transition to 
a first-order one m- The model under consideration 
here features fluctuation-induced attractive inter-vortex 
forces. Indeed, if the Josephson coupling is set to zero, 
a fractional vortex has a logarithmically divergent en¬ 
ergy while co-directed individual vortices interact log¬ 
arithmically at large distances. For a two-dimensional 


FIG. 5. Scaling behavior of the heat capacity (left column) 
and energy cumulant (right) for rj = 10“® (with A < blue 
lines) and rj = 10°'® (with A > ^j, red lines), indicates a first- 
order and a continuous phase transition respectively. Shaded 
regions around curves correspond to estimated errors. 


cross-section of a pair of fractional vortices in the two- 
component model with and phase 

winding in the first condensate, the interaction energy 

[HIES] is: 

log ^ -b TT \^|;\^Ko{r/X). (3) 

The interaction between (1,0)- and (0, l)-vortices con¬ 
tains an attractive logarithmic part and a repulsive ex¬ 
ponentially screened part: 

^("o)-s(o.i) = log ^ + TT\tlj\'^Ko{r/X). (4) 

In the absence of fluctuations, a composite vortex line 
is an axially symmetric object with finite energy per unit 
length. The attractive and repulsive forces cancel in the 
small separation limit mi Ea, however at separations 
larger than A a pair of co-directed fractional vortices can 
be mapped onto an electric dipole. Due to the long-range 
nature of the dipolar interaction, the thermal splitting 
of composite vortices into fractional vortices will lead to 
long-range attractive, short-range repulsive interaction. 
When a non-zero Josephson coupling is included, the in¬ 
teraction between fractional vortices is changed to lin¬ 
early attractive at distances larger than [33]. Then 
the splitting fluctuations of fractional vortices are largely 
confined within the range ^j. The Josephson coupling 
also changes the dipolar interactions: the phase differ¬ 
ence mode becomes massive and dipolar forces become 
short-range with the range again set by ^j. 

The results of our simulations suggest the scenario of 
a first-order phase transition originating in an attractive 
vortex interaction, i.e. in that case the superconducting 
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phase transition beyond mean-field approximation is as¬ 
sociated with proliferation of directed loops with induced 
dipolar attractive forces. These forces should tend to in¬ 
duce formation of polarized vortex clusters and phase 
separation. In the three-component case elementary vor¬ 
tices can be mapped onto Coulomb charges of different 
“colors” [5H1 . In that case the dipolar forces are stronger 
and we observe a more pronounced first-order phase tran¬ 
sition. On the other hand we found no signatures of 
a first-order transition where the Josephson length is 
smaller than the magnetic field penetration length and 
the model can be approximated by single-species repul¬ 
sively interacting directed loops. Note that a similar 
fluctuation-induced dipolar interaction between topolog¬ 
ical defects exists in gauge theories which break higher 
symmetry as well (t/(l) x 17(1) or SU(2)), where first- 
order phase transitions were also reported [131 [HI m- 
[22l [24]. Indeed for the 17(1) x 17(1) models the com¬ 
posite vortices are the lowest-energy topological excita¬ 
tions for any non-zero value of electric charge. Also in 
the SU{2) case one has composite vortices and Hopfions 
[331113] which should have attractive interaction. If elec¬ 
tric charge is decreased, the length scale at which the 
dipolar interaction sets in is increased, so in this scenario 
it requires a larger system to detect first order phase tran¬ 
sition at low electric charge. Yet, the composite vortices 
are the lowest energy topological defects at any non-zero 
value of electric charge and cannot be a priori neglected. 
The energy of topological defects is indeed different in 
models that have global symmetry, such as 17(1) x 17(1) 
which can have composite vortices due to dissipationless 
drag interaction. There, in contrast to the gauge the¬ 
ories, despite the existence of a paired phase HH H5] . 
there is indeed a tricritical point and a continuous phase 
transition at low inter component coupling [2TJ UHl |46| . 

In conclusion, we have reported that in the Lon¬ 
don limit 17(1) multiband superconductors can have a 
fluctuation-induced first-order phase transition in zero 
external field, in contrast to the inverted-XY phase 
transitions in single-band 17(1) London models [3H5|. 
We argued that the mechanism responsible for driving 
the phase transition to first-order is fluctuation-induced 
dipolar interactions between composite vortices. The 
mechanism should also apply for other multicomponent 
gauge theories including theories with higher broken sym¬ 
metry. 

Details of the numerical discretization 

The system is discretized onto a square three- 
dimensional lattice L X L X L with isotropic grid spacing 
h and periodic boundary conditions imposed in all direc¬ 
tions. Every site i has the phases and the three 

components Aiy and Ai^ of the vector (A)^. The 
phase gradients are discretized by the finite difference ap¬ 


proximation = (0-+^ — d-“^)//i with ^ = x,y, z, 

so that the kinetic energy density on site i can be written 
+ hA^y)‘^/h'^. The magnetic field energy 
density term on site i is calculated from the definition of 
a curl as an infinitesimal circulation. Using the notation 
x' = y, y' = z and z' = x, we can write the circula¬ 
tion around the plaquette P with corner in i, with area 
h and normal /i as ((V x A) • jl)i « {h~‘^ j>p A ■ dr)i = 
h {hAiyf hA^j^pf p" pf HAi^ff^. Denoting 

Aip' -f pf Aiyt! = (A X A]j^ we obtain 

A2 ((V X Af)^ = A2 E^((V X A) . A)f = [X/hf X 

hA]fp/h'^. We may furthermore rewrite the Josephson 
term as cos(6>(^) — 9^^^) = (h/^j)^ cos{9^^'> — 9^'^^)/h?. 

Since all terms in the Hamiltonian contains the factor 
l//i^, both length scales appear in the units of /i, and 
we can absorb h into A, we may set h = 1 without loss 
of generality. After rewriting the kinetic energy terms 
into cosine AY-terms (cos(a;) « I — a:^/2) we obtain the 
Hamiltonian ([^. 
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